Entanglement spectrum and boundary theories with projected entangled-pair states 
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In many physical scenarios, close relations between the bulk properties of quantum systems and 
theories associated to their boundaries have been observed. In this work, we provide an exact duality 
mapping between the bulk of a quantum spin system and its boundary using Projected Entangled 
Pair States (PEPS). This duality associates to every region a Hamiltonian on its boundary, in such 
a way that the entanglement spectrum of the bulk corresponds to the excitation spectrum of the 
boundary Hamiltonian. We study various specific models, like a deformed AKLT [T], an Ising-type 
[2], and Kitaev's toric code [3], both in finite ladders and infinite square lattices. In the latter case, 
some of those models display quantum phase transitions. We find that a gapped bulk phase with 
local order corresponds to a boundary Hamiltonian with local interactions, whereas critical behavior 
in the bulk is reflected on a diverging interaction length of the boundary Hamiltonian. Furthermore, 
topologically ordered states yield non-local Hamiltonians. As our duality also associates a boundary 
operator to any operator in the bulk, it in fact provides a full holographic framework for the study 
of quantum many-body systems via their boundary. 



PACS numbers: 

I. INTRODUCTION 

It has long been speculated that the boundary plays a 
very significant role in establishing the physical proper- 
ties of a quantum field theory. This idea has been very 
fruitful in clarifying the physics of the fractional quan- 
tum Hall effect, and is also the origin of the holographic 
principle in black hole physics. An explicit manifestation 
of this fact is the so-called area law. The area law states 
that for ground (thermal) states of lattice systems with 
short-range interactions, the entropy (quantum mutual 
information) of the reduced density operator p^, corre- 
sponding to a region A, is proportional to the surface 
of that region, rather than to the volume, at least for 
gapped systems j^HZj. Criticality may reflect itself by 
the appearance of multiplicative and/or linear logarith- 
mic corrections to the area law [8, 9 . 

Apart from the deep physical significance of this law, 
it has important implications regarding the possibility 
of simulating many-body quantum systems using tensor 
network (TN) states [T0HT3] . For instance, it has been 
shown [14] that any state of a quantum spin system ful- 
filling the area law in one spatial dimension (including 
logarithmic violations) can be efficiently represented by 
a matrix product state (MPS) [T5J[TB], the simplest ver- 
sion of a TN. 

Very recently, another remarkable discovery has been 
made with relation to the area law [17] . It has been shown 
that for certain models in two spatial dimensions, the re- 
duced density matrix of a region A has a very peculiar 
spectrum, which is called the "entanglement spectrum": 
by taking the logarithm of the eigenvalues of pA, one ob- 
tains a spectrum that resembles very much the one of a 
1-dimensional critical theory (i.e. as prescribed by con- 
formal field theory). This has been established for dif- 



ferent systems as diverse as gapped fractional quantum 
Hall states [17] or spin- 1/2 quantum magnets [18]. Inter- 
estingly, the correlation length in the bulk of the ground 
state can be naturally interpreted as a thermal length in 
one dimension [T8] . 

This is all very suggestive for the fact that the reduced 
density matrix is the thermal state of a 1-dimensional 
theory. However, there is a clear mismatch in dimen- 
sions: the Hilbert space associated to pA has two spatial 
dimensions, while the 1-dimensional theory obviously has 
only 1. Intuitively, this is clear as all relevant degrees of 
freedom of pA should be located around the boundary of 
region A. The main question addressed in this paper is 
to explicitly identify the degrees of freedom on which this 
1-dimensional Hamiltonian acts. 

We show that projected entangled-pair states (PEPS) 
[T9] give a very natural answer to that question. The de- 
grees of freedom of the 1-dimensional theory correspond 
to the virtual particles which appear in the valence bond 
description of PEPS, and that "live" at the boundary 
of region A [19] [20]. More specifically, PEPS are built 
by considering a set of virtual particles at each node of 
the lattice, which are then projected out to obtain the 
state of the physical spins. As we show, the boundary 
Hamiltonian can be thought of as acting on the virtual 
particles that live at the boundary of region A. Further- 
more, we will present evidence that, for gapped systems, 
such a boundary Hamiltonian is quasi-local (i.e. con- 
tains only short-range interactions) in terms of those (lo- 
calized) virtual particles. As a quantum phase transition 
is approached, the range of the interactions increases. 
Finally, we will show that the interactions lose their lo- 
cal character for the case of quantum systems exhibiting 
topological order. We will also show how operators in 
the bulk can be mapped to operators on the boundary. 
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The fact that the boundary Hamiltonian is quasi-local 
has important implications for the theory of PEPS which 
go well beyond those of the area law. While PEPS are 
expected to accurately represent well the low energy sec- 
tor of local Hamiltonians in arbitrary dimensions [21] , it 
has not been proven that one can use them to determine 
expectation values in an efficient and accurate way. For 
that, one has to contract a set of tensors, a task which 
could in principle require exponential time in the size 
of the lattice. In order to circumvent this problem, a 
method was introduced [19] which successively approx- 
imates the boundary of a growing region by a matrix 
product density operator, which is exactly the density 
matrix of local virtual particles discussed before. It is 
not clear a priori to which extent that density matrix 
can be approximated by a MPS; more specifically, the 
bond dimension of that MPS could in principle grow ex- 
ponentially with the size of the system if a prescribed 
accuracy is to be reached, which would lead to an expo- 
nential scaling of the computational effort. However, that 
MPS does nothing but approximate the boundary den- 
sity operator pA for different regions A. In case such an 
operator can be written as a thermal state of a quasilocal 
Hamiltonian, it immediately follows that in order to ap- 
proximate it by a MPS one just needs a bond dimension 
that scales polynomially with the lattice size [21] . and 
thus that expectation values of PEPS can be efficiently 
determined. 



II. PEPS AND BOUNDARY THEORIES 
A. Model 

We consider a PEPS, of an N v x N h spin lattice 
in two spatial dimensions. Note that one can always find 
a finite-range interaction Hamiltonian for which is a 
ground state [2 . We will assume that we have open (pe- 
riodic) boundary conditions in the horizontal (vertical) 
direction: the spins are regularly placed on a cylinder and 
the state is translationally invariant along the verti- 
cal direction [see Fig. 0]. All spins have total spin S, 
except perhaps at the boundaries where we may choose 
a different spin in order to lift degeneracies related to the 
open boundary conditions. We will be interested in the 
reduced density operator, pi, corresponding to the spins 
lying in the first £ columns; that is, when we trace all the 
spins from column £ + 1 to Nh- 

More specifically, the effective Hamiltonian, Hi, cor- 
responding to those spins, is defined through pi = 
ex.p(—Hi)/Z£, with Zg a normalization constant. We will 
be interested not only in the entanglement spectrum [17] , 
but also in the specific form of Hi and its interaction 
length, as we will define below. 

In order to simplify the notation, it is convenient 
to label the spin indices of each column with a sin- 
gle vector. We define I n = (ii, n , i2,n 5 • • • ? ^N v ,n)i where 
ik,n = S/2, -S/2 + 1, . . . , S/2 for n = 2, . . . *N h - 1 (for 



n = 1 or n = Nh we may have different spin S). Thus, 
we can write 



|*)=^c J |7 1 ,7 2 ,...,7 iVh ). 



For a PEPS we can write 



B lN h-i 

A2 ^N h -2,Aj\ 
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Here A n = (a^ n , a 2 , n , • • • , ^ Nv ^ n ), where a ki n = 
1, 2, . . . , D with D the so-called bond dimension. Each 
of the B Tl s can be expressed in terms of a single tensor, 



tr 



N v 

IP 

Lk=l 



afe,n-i,afe,; 



(3) 



where for each value of i, a, a' , A l a a , is a D x D matrix, 
7 (the indices a and f3 correspond 



with elements A 



to the virtual particles entangled along the horizontal and 
vertical directions, respectively [19]; see Fig. [T]). For the 
first (left) and last (right) column we define L 1 and R 1 
similarly in terms of the D x D matrices l l a , and r l a ,\ 
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Thus, the tensors A, and r (for which explicit expres- 
sions will be given later on) completely characterize the 
state |^), which is obtained by "tiling" them on the sur- 
face of the cylinder. The first has rank 5, whereas the 
other two have rank 4. Here we have taken all the ten- 
sors A equal, but they can be chosen to be different if 
the appropriate symmetries are not present. 



B. Boundary density operator 

We now want to express the reduced density operator 
pi in terms of the original tensors. In order to do that, 
we block all the spins that are in the first £ columns, and 
those in the last N^ — £, and define 

L Ia = L h B l2 ...B h , R h = B Ie+1 . . . B lN h-^R lN h , 

(6) 

where we have collected all the indices 1\ , . . . , lg in I a 
and the rest in 1^. With this notation, the state can 
be considered as a two- leg ladder, ie N^ = 2, and 1 = 1, 
where pg is the density operator corresponding to a single 
leg. Thus, we have 



(7) 
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FIG. 1: Top: We consider an N v xNh spin lattice in a cylindri- 
cal geometry. The PEPS is obtained by replacing each lattice 
site with a tensor A, and contracting the virtual indices a 
and f3 along the horizontal and vertical directions. Bottom: 
We cut the lattice into two pieces, left and right. The virtual 
indices a of the tensors A along the cut are shown. The state 
acts on the spins (ik,n) as well as on the virtual spins 
along the cut. 



In terms of those operators, it is very simple to show that 

P/ = £lXr><xH (T\^la R ^l\T') (11) 
r,r' 

where |T) is an orthonormal basis of the range of <tl, g\ 
is the transpose of ol in the basis |A), and where we have 
defined an orthonormal set (in the spin space) 

Now, denning an isometric operator that transforms the 
virtual onto the spin space U = J2 F |xr)(r|, we have 




The isometry U can also be used to map any operator 
acting on the bulk onto the virtual spin space; note that 
this map is an isometry and hence not injective, i.e. a 
boundary operators might correspond to many different 
bulk operators. This is of course a necessity, as U is 
responsible for mapping a 2-dimensional theory to a 1- 
dimensional one. 



C. Boundary Hamiltonian 



It is convenient to consider the space where the vectors 
L 1 and R 1 act as a Hilbert space, and use the bra/ket 
notation there as well. That space, that we call virtual 
space, is the one corresponding to the ancillas that build 
the PEPS in the valence bond construction [15] . They 
are associated to the boundary between the £-th and the 
£ + 1st columns of the original spins. The dimension is 
thus D Nv (see Fig. [l). In order to avoid confusion with 
the space of the spins, we have used \v) to denote vectors 
on that space. We can define the (unnormalized) joint 
state for the first t columns and the virtual space, |^l) 5 
and similarly for the last columns, \^r), as 



-) = £|L'")|I a ), \* R ) = J2\R Ib )\h) (8) 



with 



il'«) = y^ia), i#<o = y^£|A), 



(9) 



and | A) the canonical orthonormal basis in the corre- 
sponding virtual spaces. The state can then be 
straightforwardly denned in terms of those two states. 
The corresponding reduced density operators for both 
virtual spaces are 



h\ 



(10) 



The previous equation shows that p£ is directly related 
to the density operators corresponding to the virtual 
space of the ancillary spins that build the PEPS. In par- 
ticular, if we have a\ = or =: (eg., when we have the 
appropriate symmetries as in the specific cases analyzed 
below), then pg = Ua^W . The reduced density opera- 
tor p£ is thus directly related to that of the virtual spins 
along the boundary. Since U is isometric it conserves the 
spectrum and thus the entanglement spectrum of pi will 
coincide with that of g\. By writing o\ = exp(— H^), we 
obtain an effective one-dimensional Hamiltonian for the 
virtual spins at the boundary of the two regions whose 
spectrum coincides with the entanglement spectrum of 
Pi- 

We will be interested to see to what extent is a 
local Hamiltonian for the boundary (virtual) space. We 
can always write as a sum of terms where different 
spin operators. For instance, for D = 2, we can take 
the Pauli operators a a (a = x, y, z) acting on different 
spins, and the identity operator on the rest. We group 
those terms into sums h ni where each h n contains all 
terms with interaction range n, i.e., for which the longest 
contiguous block of identity operators has length N v — 
n. For instance, ho contains only one term, which is 
a constant; hi contains all terms where only one Pauli 
operator appears; and h^ v contains all terms where no 
identity operator appears. We define 



la 



h 



d n = tr(^)/2^ , 



(14) 
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which expresses the strength of all the terms in the 
Hamiltonian with interaction length equal to n. A fast 
decrease of d n with n indicates that the effective Hamil- 
tonian describing the virtual boundary is quasi-local. In 
the examples we examine below this is the case as long 
as we do not have a quantum phase transition. In such a 
case, the length of the effective Hamiltonian interaction 
increases. 



D. Implications for PEPS 



of the chain, when cut into two parts, has zero dimen- 
sions, it will help us to understand the 2D systems. We 
take N v = 1 so that the PEPS reduces to a MPS. We can 
use the theory of MPS [I5j [16] to analyze the properties 
of the completely positive map (CPM) £ (the matrices A 1 
of the MPS are the Kraus operators of the CPM). In the 
limit Nh — » oo, cr& is nothing but the fixed point of such 
a CPM. For gapped systems, £ has a unique fixed point, 
and thus is unique. For gapless systems, £ becomes 
block diagonal (and thus there are several fixed points), 
the correlation length diverges, and we can write 



In case can be written in terms of a local bound- 
ary Hamiltonian one can draw important consequences 
for the theory of PEPS. In particular, it implies that the 
PEPS can be efficiently contracted, and correlation func- 
tions can be efficiently determined. The reason can be 
understood as follows. Let us consider again the cylin- 
drical geometry (Figjl]), and let us assume that we want 
to determine any correlation function along the vertical 
direction, eg at the lattice points 1) and (£,x). It is 
very easy to show that such a quantity can be expressed 
in terms of gl and or. If we are able to write these two 
operators as Matrix Product Operators (MPO), ie as 



(16) 



D 



tr[M^...M^^] \ii,...,i Nv )(ji 



,JN V 



in ijn i — 1 



(15) 

where the M' are D' x D' matrices, then the correlation 
function can be determined with an effort that scales 
as N v (D f ) 6 . It has been shown by Hastings [21] that 
if an operator can be written as exp(— iJ^/2), where 
is quasilocal, then it can be efficiently represented by an 
MPO; that is, the bond dimension D' only scales polyno- 
mially with N v . Thus, we have that the time required to 
determine correlation functions only scales polynomially 
with N v . 

Later on, when we examine various examples, we will 
use MPO to represent 07,. In that case, we can directly 
check if we obtain a good approximation by using a MPO 
just by simply observing how much errors increase when 
we decrease the bond dimension D' . We will see that 
the error increases when we approach a quantum phase 
transition. Furthermore, whenever 0-5 can be well ap- 
proximated by a MPO, we can use the knowledge gained 
in the context of MPS [15J[T6] to observe the appearance 
of a quantum phase transition in the original PEPS. For 
that, we just have to recall that the correlation length, 
£, is related to the two largest (in magnitude) eigenval- 
ues, Ai 5 2, of the matrix J2i M hl ; £ = 1/ log(|Ai/A 2 |). For 
l^i I = |^2 1, the correlation length diverges indicating the 
presence of a quantum phase transition. 



E. Qualitative discussion 

In order to better understand the structure of cr^, let us 
first consider a ID spin chain. Even though the boundary 



where B is the number blocks which coincides with the 
degeneracy of the eigenvalue of £ corresponding to the 
maximum magnitude. In such case, the weights p n de- 
pend on the tensors I and r which are chosen at the 
boundaries. For critical systems, one typically finds that 
D increases as a polynomial in N v such that one obtains 
logarithmic corrections to the area law [9] [22] . 

The 2D geometry considered here reduces to the ID 
case if we take the limit Nh — >• 00 by keeping N v finite. 
According to the discussion above, we expect to have 
a unique 0-5 if we deal with a gapped system. As we 
will illustrate below with some specific examples, this 
operator can be written in terms of a local Hamiltonian 
Hb of the boundary virtual space which is quasilocal. As 
we approach a phase transition, the gap closes and the 
correlation length diverges. In some cases, the boundary 
density operator can be written as a direct sum (16), 



eventually leading to the loss of locality in the boundary 
Hamiltonian. 



III. NUMERICAL METHODS 

In order to determine o\> we make heavy use of the 
fact that is a PEPS. We have followed three differ- 
ent complementary numerical approaches that we briefly 
describe here. 



A. Iterative procedure 

First of all, for sufficiently small values of N v (typically 
N v < 12) we can perform exact numerical calculations 
and determine (Jl,r according to (10). The main idea 
is to start from the left and find first ctl for £ — 1 by 
contracting the tensors l % appropriately. Then, we can 
proceed for i — 2 by contracting the tensors A 1 corre- 
sponding to the second column. In this vain, and as long 
as N v is sufficiently small we can determine ctl,r for all 
values of £ and Nh- 
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B. Exact contractions and finite size scaling 

The second (exact) method is a variant applicable to 
larger values of N v (typically up to N v = 20) but re- 
stricted to a finite width in the horizontal direction. It 
consists in exactly contracting the internal indices of two 
adjacent blocks of size N v /2 x N^. These two blocks are 
then contracted together in a second step. Although lim- 
ited by the size 2 Nv + 2Nh of the half-block (which has to 
fit in the computer RAM) , this approach can still handle 
systems of size 20 x 2 or 16 x 8 and be supplemented by 
a finite size scaling analysis. 



C. Truncation method 



Finally, to take the Nh — >• oo limit we can use the meth- 
ods introduced in [19] to approximate the column oper- 
ators. The main idea is to represent those operators by 
tensor networks with the structure of a MPS. We contract 
one column after each other, finding the optimal MPS 
after each contraction variationally. In particular, since 
we will consider translationally invariant states, we can 
choose the matrices of the corresponding MPS all equal, 
which simplifies the procedure. We can even approach 
the limit N v ,Nh oo as follows (see also [I2j|23]): (i) 
we start out with i = 1, and contract the second column, 
obtaining another tensor network with the same MPS 
structure, but with increased bond dimensions, (ii) We 
continue adding columns, up to some £ = r, where we 
start running out of resources. At that point, we have a 
tensor network with the MPS structure representing ctl- 
Let us denote by o the basic tensor of that network, 
where n = 1, . . . , D 2 and a, /3 = 1, . . . , D 2r (n denotes the 
index in the horizontal direction), (iii) When the bond 
indices a, j3 grow larger than some predetermined value, 
say D c < D 2r we start approximating the tensor network 
by one with bond dimension D c as follows. We first con- 
struct the tensor K a ^.^^> = "}2 n C™^C™,^,. Later on 
we will always deal with the case in which K is hermitean 
(when considered as a matrix); if this is not the case, one 
can always choose a gauge where it is symmetric [IB] . We 
determine the eigenvector, Xp^p, corresponding to the 
maximum eigenvalue of diagonalize X, consider the 
D c largest eigenvalues and build a projector onto the cor- 
responding eigenspace. We then truncate the indices a 
and j3 by projecting onto that subspace. (iv) We continue 
in the same vein until the truncated tensor structure con- 
verges, which corresponds to the limit — > oo. (v) We 
can do the same with cjr by going from right to left. For 
the examples studied below, gl = <r# =: = ', and 
thus we just have to carry out this procedure once. 



IV. NUMERICAL RESULTS FOR AKLT 
MODELS 

We now investigate some particular cases. We concen- 
trate on the AKLT model [I] |24], whose ground state, 
can be exactly described by a PEPS with bond di- 
mension D = 2, as shown in Figs. [2] and and [3] The spin 
in the first and last column have S = 3/2, whereas the 
rest have S — 2. The AKLT Hamiltonian is given by a 
sum of projectors onto the subspace of maximum total 
spin across each nearest neighbor pair of spins, 

#AKLT= P n%> ( 17 ) 

<n,m> 

(s) 

where PA,m is the projector onto the symmetric subspace 
of spins n and m. This Hamiltonian is su(2) and transla- 
tionally invariant. This invariance is inherited by the vir- 
tual ancillas, and thus and will also be. These sym- 
metries can be used in the numerical procedures. Note 
that if Hi) has this symmetries and has short-range inter- 
actions, then since the ancillas have spin 1/2 (as D = 2), 
it will be generically critical. 




-4- -K-k 

FIG. 2: (Color online) (a) Ribbon made of two (N h = 2) 
coupled periodic S=l/2 Heisenberg chains (2-leg ladder), (b) 
Groundstate of a 2-leg S=3/2 AKLT ladder. Each site is split 
into three spins- 1/2 (red dots). Nearest neighbor spins- 1/2 
are paired up into singlet valence bonds, (c) PEPS represen- 
tation for S=3/2 and S=2 sites of AKLT wavefunctions in the 
valence bond (singlet) picture (for connection to the "maxi- 
mally entangled picture" see text). Open squares stand for 
the r™ lja2)OL3 and A^ ljQ , 2)Q;3jQ!4 tensors defined in the text and 
open circles correspond the to 2 x 2 matrix [0, 1; —1, 0]. 

The lattice is bipartite. It is convenient to apply the 
operator ex.p(i7rS y /2) to every spin on the B sublat- 
tice: this unitary operator does not change the prop- 
erties of pi but slightly simplifies the description of the 
PEPS. Thus, we can write the AKLT Hamiltonian as 
in (|TtJ) but with P n?m := exp(z7r(/ n 6' 2y , n + 

ex.p(i7r(f n Sy }n + fmSy ym ), with f n = 0, 1/2 



6 



if the spin n is in the A or B sublattice, resp. 

We will study finite A^-legs ladders, as well as infinite 
square lattices. We will start out in the next subsection 
with the simplest case of Nh = 2. Note that for this 
particular case the subsystem we consider when we trace 
one of the legs is a spin chain itself, so that density oper- 
ator p£ = i already describes a 1-dimensional system and 
thus the physical spins already represent the boundary. 
In such a case, we do not need to resort to the PEPS for- 
malism but we can also study other model Hamiltonians 
besides the AKLT one. For example, we will consider 
the 5ii(2)-symmetric Heisenberg ladder Hamiltonian of 
S = 1/2 [Fig. §a)] 



Heis 



= £ 



(18) 




(a) 



N h /2 



b— 



d — 



(b) 



where the exchange couplings J ny7n are parametrized by 
some angle #, i.e. J\ eg = cos 6 (J run g = sin#) for nearest- 
neighbor sites n and m on the legs (rungs) of the ladder. 
Although the ground state has no simple PEPS repre- 
sentation, it can be obtained numerically by standard 
Lanczos exact diagonalization techniques on finite clus- 
ters of up to 14 x 2 sites [T8] . Similarly to the AKLT 
2-leg ladder [Fig. [2^b)] , it possesses a finite magnetic 
correlation length £ which diverges when 6 — > (decou- 
pled chain limit). The opposite limit = tt/2 (0 = —tt/2) 
corresponds to decoupled singlet (triplet) rungs (strictly 
speaking, with zero correlation length). 

For infinite systems, we will also be interested in the 
behavior of along a quantum phase transition. To 
this aim, we will also consider a distorted version of the 
AKLT model, and define a family of Hamiltonians 

H(A)= Q„(A)Q m (A) J P„, m Q„(A)Q m (A), (19) 



where Q n (A) = e~ 8ASz > n . Note that the Hamiltonian is 
translationally invariant and has u(l) symmetry. As A 
increases, it penalizes (nematic) states with S z = 0, and 
thus the spins tend to take their maximum value of S%. 
As we will show, there exists a critical value of A where 
a quantum phase transition occurs. 



A. 2— leg ladders : comparison between AKLT and 
Heisenberg models 

Let us start out with the 5^(2)-symmetric A = 
AKLT model in a two-leg ladder configuration, where 
pi corresponds to state of one of the legs; that is, we 
take Nh = 2, t = 1, and all spins have S = 3/2 as 
shown in Fig. [5|b). The Hamiltonian is gapped [TJ [24] . 
and the ground state is a PEPS with bond dimension 
D = 2. The tensors corresponding to the two legs, I and 
r, coincide and are given by r™ >Q!2)Q!3 = (s rn \a 1: a 2 ,a s }, 
where = ±1/2, and \s m ) is the state in the symmet- 
ric subspace of the three spin 1/2 with S z \s m ) = m|s m ), 
m = -3/2,-1/2,1/2,3/2. 



FIG. 3: (Color online) (a) 4-leg (N h = 4) AKLT ladder on a 
cylinder partitioned (dotted green line) into two halves, (b) 
Schematic representation of the density matrix a\ of a 4-leg 
(Nh = 4, £ = 2) AKLT ladder. After being "cut" the two 
halves are "glued" together (physical indices are contracted). 



We first examine the entanglement spectrum of 
computed on a 16 x 2 ladder. It is shown on Fig.Qb) as a 
function of the momentum along the legs, making use of 
translation symmetry (the vertical direction is periodic) 
enabling to block-diagonalize the reduced density matrix 
in each momentum sector K. Note that it is also easy 
to implement the conservation of the z-component S z of 
the total spin so that each eigenstate can also be labelled 
according to its total spin S. The low-energy part of the 
spectrum clearly reveals zero-energy modes at K = and 
K = 7r consistent with conformal field theory of central 
charge c = 1. 

It is of interest to compare the 2-leg AKLT results 
to the ones of the 2-leg S=l/2 Heisenberg ladder jT8[ ) 
sketched in Fig.[2|a) and investigated in Ref.[18j Fig.Qa) 
obtained on a 14 x 2 ladder for a typical parameter 
= 7r/3 shows the entanglement spectrum of p£ which, 
again, is very similar to that of a single nearest-neighbor 
Heisenberg chain. As mentioned in Ref. [18j in first ap- 
proximation, varying the parameter (and hence the lad- 
der spin-correlation length) only changes the overall scale 
of the energy spectrum. Hence, it has been suggested [18] 
to connect this characteristic energy scale to an effective 
inverse temperature /3 e ff. 

The above results strongly suggest that is " close" to 
a one-dimensional nearest-neighbor Heisenberg Hamilto- 
nian. To refine this statement and make it more precise, 
we perform an expansion in terms of sii(2)-symmetric 
extended-range exchange interactions, 



H h = AM, 



• ^ A r Sk • S/e +r + RX . 



(20) 



where RX stands for the "rest", i.e. (small) multi-spin 
interactions. The amplitudes A r can be computed from 
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FIG. 4: (Color on line) Entanglement spectra of if 6 (w.r.t. the 
groundstate energy £ ) versus total momenta K in the chain 
(vertical) direction, (a) 2-leg (14 x 2) quantum Heisenberg 
ladder, (b) 2-leg (16 x 2) AKLT ladder and (c) 8-leg (16 x 8) 
AKLT ladder. The eigenvalues are labelled according to their 
total spin quantum number using different symbols (according 
to the legend on the graph). 



simple trace formulas, 



A r = —tr{H b J2* z k Vk +r }/2 N ' J 



(21) 



requiring the full knowledge of the eigenvectors of (i.e. 
of <7b). Ao is fixed by some normalization condition, e.g. 
tr<j5 = 1. Assuming X is normalized as an extensive 
operator in N v , i.e. ^j-tr{X 2 } = 2^, the amplitude R is 
given by: 



R* = ^tr{^ 2 }/2^ - N v Al - ± A " ' ( 22 ) 

v r—l 



The coefficients A r and R of 2-leg Heisenberg ladders 
are plotted in Fig. [5|a) as a function of the parameter 

0, both in the Haldane (J rU ng < i.e. ferromagnetic) 
and rung singlet phases (J run g > i.e. antiferromag- 
netic). Generically, we find that is not frustrated, 

1. e. all couplings at odd (even) distances are antiferro- 
magnetic (ferromagnetic), A r > (A r < 0). Clearly, 
the largest coupling is the nearest-neighbor one (r = 1). 
Fig. [5|b) shows the relative magnitudes of the couplings 
at distance r > 1 w.r.t. A\. These data suggest that 
the effective boundary Hamiltonian is short range, 
especially in the strong rung coupling limit (0 — >• 7r/2) 
where \A r > j A r \ — >• for r' > r. The amplitude A\ of 
the nearest-neighbor interaction can be identified to the 
effective inverse temperature f3 e fi which, therefore, van- 
ishes (diverges) in the strong (vanishing) rung coupling 
limit. 



FIG. 5: (Color on line) (a) Amplitudes A r of the (isotropic) 
spin-spin couplings up to distance r = 7 of the effective 
boundary Hamiltonian of a quantum Heisenberg 2-leg lad- 
der in the Haldane and rung singlet phases vs 0. (b) Ratio 
of the same amplitudes normalized to the nearest-neighbor 
coupling (r = 1). Computations are carried out on 12 x 2 
(open symbols) and 14 x 2 (closed symbols) systems. Note 
that when — » tt/2 (decoupled rung singlets), A r /Ao — > 
for r > 1 and all the weights of the reduced density matrix 
becomes equal to 2~ Nv (Ao — In 2). 
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FIG. 6: (Color on line) 2-leg quantum Heisenberg ladders - 
Ratio of the amplitudes \A r \ by the nearest-neighbor ampli- 
tude A\ plotted using a logarithmic scale as a function of r for 
different values of 0. (a) Antiferromagnetic and (b) ferromag- 
netic leg couplings (the rung couplings are antiferromagnetic 
in both cases). 



Next, we investigate the functional form of the decay 
of the amplitudes \A r \ with distance. The ratio |A r |/Ai 
versus r are plotted (using semi-log scales) in Figs.[6ja,b) 
for 12 x 2 and 14 x 2 Heisenberg ladders with differ- 
ent values of 0. Similar data for a 20 x 2 AKLT lad- 
der is shown in Fig. (7^a), providing clear evidence of 
exponential decay of the amplitudes with distance, i.e. 
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\A r \ ~ exp (— r/£b). The Heisenberg ladder data are also 
consistent with such a behavior (even though finite size 
corrections are stronger than for the AKLT case, espe- 
cially when — >• or n). It is not clear however how deep 
the connection between the emerging length scale £5 and 
the 2-leg ladder spin correlation length £ is. Note that 
the latter can be related [18 to some effective thermal 
length associated to the inverse temperature /3 e ff oc A\. 

Thanks to the PEPS representation of their ground 
state, AKLT ladders can be (exactly) handled up to 
larger sizes than their Heisenberg counterparts (typically 
up to N v = 20) enabling a careful finite size scaling anal- 
ysis of the boundary Hamiltonian (20). As shown in 
Fig. [8ja), we observe a very fast (exponential) conver- 
gence of the coefficients A r with the ladder length N v . 
Hence, one gets at least 7 (3) digits of accuracy for all 
distances up to r = 5 (r = 7). 
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B. N h -\eg AKLT ladder 

Now we consider the AKLT model on an iV^-leg lad- 
der configuration; we take i — Nh/2. The spins in the 
first and last legs have S = 3/2, and the correspond- 
ing tensors coincide with the ones given above. The rest 
of the spins have 5 = 2, and the corresponding tensor 
is ^ai,a 2 ,a 3 ,a4 = ( s m |«i , <*2 , <*3 , a 4 ) , where a* = ±1/2, 
and \s m ) is the state in the symmetric subspace of the 
four spin 1/2 with S z \s m ) = m|s m ), m = —2,-1,0,1,2 
(see Fig. [2jc)). An example of a 4-leg ladder and of a 
schematic representation of p£ is shown in Fig. [3j 
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FIG. 8: (Color on line) (a) Finite size scaling of the ampli- 
tudes A r for a 2-leg AKLT ladder vs 1/N V (N v = 14, 16, 18, 
20). (b) Finite size scaling of the amplitudes A r for iV^-leg 
AKLT ladders vs 1/Nh at fixed N v = 16 (open symbols) or 
N v = 18 (filled and + symbols), (c) VN entropy per unit 
length (normalized by ln(2)) vs 1/Nh at fixed N v = 16. 



FIG. 7: (Color on line) AKLT ladders - (a) Ratio \A r \/At 
plotted using a logarithmic scale as a function of r. Results 
are approximation-free for finite Nh while the Nh oo limit 
is obtained by finite size scaling (see Fig.|8jb)). (b) Compar- 
ison with y/d r +i/d2 (full symbols) computed (see text) on 
2-leg and infinitely long (Nh = oo) cylinder. 



In fact, as pointed out previously, the boundary Hamil- 
tonian Hh should not contain only two-body spin inter- 
actions. However, the total magnitude of all left-over 
(multi-body) contributions, i?, is remarkably small in the 
AKLT 2-leg ladder : as shown in Fig. (8^a), R < A 4 . In 
fact, the full magnitude of all man y-body terms extend- 
ing onr + 1 sites is given by ^^r+i and can be compared 
directly to \A r \ (after proper normalization). Fig. [7jb) 
shows that and |A r |/Ai are quite close, even 

at large distance. Note however that multi-body interac- 
tions are significantly larger in the boundary Hamiltonian 
of the Heisenberg ladder, as shown in Fig [5] (although no 
accurate finite size scaling analysis can be done in that 
case). 



Let us now follow the same analysis ( 20 ) of the bound- 
ary Hamiltonian as we did for the case of 2 legs. The 
decay with distance of the coefficients A r are reported 
in Fig. [7ja) for 4-leg, 6-leg and 8-leg AKLT ladders. 
Clearly, the decay is still exponential with distance for all 
values of Nh studied but the characteristic length scale 
associated to this decay (directly given by the inverse of 
the slope of the curve in such a semi-log plot) smoothly 
increases with Nh- A careful finite size scaling is per- 
formed in Fig. |8jb) to extract the Nh —> oo limit of 
all A r (accurate up to r = 7). The extrapolated val- 
ues are reported in Fig.[7ja) showing that A r also decays 
exponentially fast with r in an infinitely long cylinder 
(Nh = oo). The characteristic emerging length scale is 
estimated to be still very short around 1. 

Lastly, we compute the Von Neumann entanglement 
entropy defined by Sy-^(pi) = ~^{pi m Pi} with the nor- 
malization tr pi = 1. 5Vn scales like N v ("area" law) and 
is bounded by N v In 2. Fig. |8jc) shows that the entropy 
converges very quickly with Nh to its thermodynamic 
value which is very close to the maximum value. The 
entanglement of the two halves of the AKLT cylinder is 
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therefore very strong. 



C. Thermodynamic limit and phase transitions 

Now we consider the N v ,Nh — ^ oo for the deformed 
AKLT model in order to investigate the phase transition. 
We will compare some of the results with the 2-leg ladder 
as well. The spins in the first and last legs have S = 
3/2, and the rest S — 2. The corresponding tensors are 
defined according to 



r 



(s m \Q(A)\a 1 ,a 2 ,a 3 ), 



A: 



(5 m |Q(A)|ai,a2,a 3 ,a 4 ), 



(23) 



where = ±1/2, and \s m ) is the state in the symmet- 
ric subspace of the three (four) spin 1/2 with S z \s m ) = 
m\s m ), m = -3/2,-1/2,1/2,3/2 (m = -2,-1,0,1,2), 
respectively. 

We will use the approximate procedure sketched in Sec- 
tion III-C. In particular, for N v larger than the correla- 
tion length the obtained tensors C™ o will be independent 
of N v . We have considered those tensors (with D c = 50 
and 100 iterations), and built and out of them. 
Note that the su(2) symmetry is explicitly broken by a 
finite A so that it becomes more convenient to use the 
variable d n of Eq. (14) instead of A r to probe the spatial 
extent of H^. We recall that (dn) 1 / 2 is the mean ampli- 
tude of all interactions acting at distance r = n — 1. We 
have plotted in Fig. [9] all d n , n < N v /2, for N v = 16 as 
a function of A. As A increases, we see that the inter- 
action length of the effective Hamiltonian increases and 
one sees a long-range interaction appearing. This indi- 
cates that we approach a phase transition. For the case 
of the ladder, the interaction length remain practically 
constant for the same range of variation of A. 

Similarly to the investigation of the Heisenberg lad- 
der [T8] . it is interesting to define an effective inverse 
temperature via the amplitude of the nearest-neighbor 
interaction, 



0, 



eff 




(24) 



where the pre-factor is introduced conveniently so that 
/3 e ff = A\ in the sw(2)-symmetric limit A = 0. As seen 
in the inset of Fig. [9j the inverse temperature of the lad- 
der scales linearly with A. For the infinite cylinder, no 
singularity of /3 e ff is seen at the cross-over between short 
and long-range interactions. 

Next, we plot the inverse correlation length as a func- 
tion of A both for one dimension (i.e. an infinitely long 
ladder) and for two dimensions (i.e. N v = Nh = oo) in 
Fig. [ToFa), obtained with D c = 150 and 100 iterations 
(no difference are observed by taking D c = 50 and 50 
iterations). Clearly, the divergence of £ (i.e. £ _1 =0) 
shows the appearance of a phase transition at A = 0.0061 
in two dimensions. In contrast, £ _1 never crosses zero 
in the case of the ladder (i.e. in one dimension). We 




FIG. 9: (Color on line) AKLT model with finite "nematic" 
field A - (a) Relative amplitude y/d r +i/d2 in a 2-leg ladder 
plotted using a logarithmic scale as a function of r; (b) same 
for an infinitely long cylinder (Nh = oo). From bottom to 
top, A is incremented from to A max by constant steps, (c) 
Inset : effective temperature /3 e ff (see Eq. 24 ) versus A for the 
two cases reported in (a) and (b). All results are obtained for 
N v = 16 (D c = 50, and 100 iterations such that the tensors 
C already converge). 



have compared £ with the "emerging" length scale £5 ob- 
tained by fitting the decay of the coefficients of as 
yjd r +\jd<i ~ exp (— r/£b) on N v = 16 2-leg and infinitely 
long (i.e. Nh = 00) cylinders. In the two leg ladder, we 
see that the divergence of the correlation length £ for 
A — »> 00 results from the interplay between (i) a (moder- 
ate) increase of the range £5 of the Hamiltonian and 
(ii) a linear increase with A of the effective temperature 
scale /3 e ff, therefore approaching the T e ^ — >> limit when 
A 00. This contrasts with the case of two dimensions 
(N v = Nh = 00) where the divergence of £ occurs at fi- 
nite effective temperature when becomes "sufficiently" 
long-range. It is however hazardous to fit the decay of 
the coefficients of to obtain its functional form at the 



phase transition. Finally, in Fig. 101)) we have plotted 
the truncation error made by taking different D c in the 
limit Nh — >• 00, and, again around A w 0.006 the error 
increases. This is consistent with the expectation that as 
Hi) contains longer range interaction, the boundary den- 
sity operator 0-5 requires a higher bond dimension to be 
described as a TN state. 



V. NUMERICAL RESULTS FOR ISING PEPS 

We now continue by considering the Ising PEPS intro- 
duced in [2]. They all have bond dimension D = 2 and 
exhibite the Z2-symmetry of the transverse Ising chain. 
They depend on a single parameter, G [0,7r/4]. For 
9 ~ 7r/4 one has a state with all the spins pointing in 
the x direction, whereas for 9 ~ the state is of the 
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FIG. 10: (Color online) (a) Inverse correlation length £ _1 
vs A for both 2-leg ladder and infinitely long cylinder (Nh = 
oo). These data correspond to the infinite circumference limit, 
i.e. N v = oo. The arrow marks the phase transition in the 
infinitely long cylinder. Comparison with the inverse of the 
"emerging" length scale £b obtained bv fitting the decay of 
the coefficients of Hb plotted in Fig. k)[ a) as y/d r +i/d,2 ~ 
exp (— r/£b)- (b) Truncation error in the Nh —> oo procedure. 
The results are compared with those obtained with D c = 150, 
and 100 iterations have always been used. 

GHZ type (a superposition of all spins up and all down). 
In the thermodynamic limit (N v ,Nh — » oo) for ~ 0.35 
they feature a phase transition, displaying critical behav- 
ior, where the correlation functions decay as a power law. 
Thus, by changing we can investigate how the bound- 
ary Hamiltonian behaves as one approaches the critical 
point. 



A. 2— leg ladders 

The tensors corresponding to the two legs, I 
and r, coincide and are given by r™ a2 a = 
^m{ a i) a m{ a 2)^m{ a ?>)i where m = 0, 1, a t = ±1/2 and 
a m (a) are parametrized as an(— 1/2) = ai(l/2) = cos# 
and a (l/2) = a^-1/2) = sin (9. 

As seen in Fig. [TT] the entanglement spectrum of the 2- 
leg ladder is gapped for all values and resemble the one 
of an Ising chain (equally spaced levels) with small quan- 
tum fluctuations revealed by the small dispersion of the 
bands. The effective inverse temperature, qualitatively 
given by the gap (or the spacing between the bands), 
decreases for increasing 0. 

The interaction length of the boundary Hamiltonian 
for the ladder is displayed in Fig. [l2^a). The strength 
of the interactions decay exponentially with the distance 
for all values of 6. As we increase this angle, one only 
observes a decrease of the interaction length. Note that 
as opposed to the AKLT models studied in the previous 
sections, d± ^ 0. Indeed, there always exists a term 



FIG. 11: Entanglement spectrum of a 16 x 2 Ising PEPS lad- 
der versus momentum along the ladder leg direction. Com- 
parison between — 0.2 (a) and — 0.5 (b) using the same 
energy scale, (c) Zoom in of the low-energy part of (b). 

with a single Pauli operator a X: describing an effective 
transverse field in H^. Thus, that Hamiltonian is given 
by a transverse Ising chain in the non-critical region of 
parameters. 

We have also plotted the inverse correlation length £ -1 
as a function of in Fig. [l3] (blue empty dots). While 
the correlation length increases as decreases, it only 
tends to infinite in the limit — » 0, as it must be for a 
GHZ state. No signature of a phase transition is found 
otherwise. 



B. Thermodynamic limit and phase transitions 

We now move to the case of an infinitely long cylinder. 
As above to grow the cylinder in the horizontal direction, 
one considers rank-5 tensors, which here take the form 
^ai,a 2) a 3 ,a 4 = a ™ )a m (a 2 )a m (a 3 )a m (a 4 ) , and use the 
same approximation scheme with 100 iterations as before. 

The parameters d n describing the boundary Hamilto- 
nian Hi) behave very differently in the ladder and infinite 
cylinders as shown in Fig. [12] While for the Ising PEPS 
ladder H b remains short-ranged with exponential decay 
of d n vs n, the infinite cylinder shows a transition to- 
wards long-range interactions suggesting the existence of 
a phase transition. This is very similar to what occurred 
in the AKLT distorted model. 

As long as H b remains short-range, the density matrix 
Pi can be (qualitatively) mapped onto the thermal den- 
sity matrix of an effective quantum Ising chain (including 
a "family" of transverse-like fields) and, therefore, no or- 
dering is expected (at finite effective temperature). A 
phase transition however can appear when H b becomes 
long-ranged as it is the case for an infinitely long cylin- 
der. This is evidenced by the the behavior of correla- 
tion lengths computed for the 2-leg and infinitely long 
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FIG. 12: (Color online) Ising PEPS - Relative amplitude 
y/dn/di in a 2-leg ladder (a) and in an infinitely long (Nh = 
oo) cylinder (b) as a function of n for varying from 0.1 to 
- 0.51 with constant intervals (0.1, 0.1585, 0.2171, 0.2756, 
0.3342, 0.3927, 0.4512 and 0.5098, from top to bottom). Log- 
arithmic scales are used on the vertical axis in both (a) and 
(b). Inset: Ratio of the effective transverse field y/di over 
the effective Ising nearest-neighbor coupling \fd^ versus 6 for 
Nh = oo (D c = 50 and 100 iterations). All results are ob- 
tained for N v — 12. 



(Nh = oo) cylinder and reported in Fig. 13 These corre- 
lation lengths are compared to the respective " emerging" 
length scales characterizing the decay of V^n with n. 
In the 2-leg ladder case, £5 increases quite moderately 
when — )> (£& ~ 1) so that the divergence of the corre- 
lation length £ in this limit is only attributed to a van- 
ishing of the effective temperature scale T e ^ . In contrast, 
as for the AKLT model, the phase transition in two di- 
mensions occurs at finite (effective) temperature at the 
point where £5 — >• 00. 

In summary, these results evidence that whenever we 
approach a phase transition, the interaction length of the 
boundary Hamiltonian increases. 



FIG. 13: (Color on line) Ising PEPS - Inverse correlation 
length £ _1 vs A for both 2-leg ladder and infinitely long cylin- 
der (Nh = 00, D c = 150 and 100 iterations). These data cor- 
respond to the infinite circumference limit, i.e. N v = 00. The 
arrow marks the phase transition in the infinitely long cylin- 
der. Comparison with the inverse of the "emerging" length 
scale £ b " 1 obtained by fitting the decay of the coefficients plot- 
ted in Fig. 12 as \fa\ ~ exp (- 
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the spins in each block by one tensor of the form 



/pi,2,*2,3,*3,4,U,: 
^ 1 Q!l,Q;2,Q!3 5 Q! 4 



if i x , x +i — a x+i — a x mod 2 
otherwise. 

(26) 

Here, i x ,x+i denotes the spin located between the 
bonds a x and a x +\ (numbered clockwise) as shown in 
Fig, p^b). It can be checked straightforwardly that the 
resulting tensor network is an eigenstate of the Hamilto- 
nians of Eq. (25). Excitations of the model correspond 
to violations of fox-terms (charges) or fo^-terms (fluxes), 
which always come in pairs [3]. 

We put the code state on a cylinder of Nh x N v ten- 
sors (i.e., 2Nh x 2N V sites), where we choose boundary 
conditions 



VI. TOPOLOGICAL KITAEV CODE 

Let us finally consider systems with topological order. 
We will focus on Kitaev's code state [3]: It can be defined 
on a square lattice with spin-^ systems (qubits) on the 
vertices, with two types of terms in the Hamiltonian, 



X 



04 



z m 



(25) 



(where X and Z are Pauli matrices), each of which acts 
on the four spins adjacent to a plaquette, and where the 
hx and hz form a checkerboard pattern (see Fig.[l4^a)). 
The ground state subspace of the code state can be rep- 
resented by a PEPS with D = 2 [2]; a particularly con- 
venient representation is obtained by taking 2x2 blocks 
of spins across hz type plaquettes, and jointly describing 



|x*)=cosf |0)« 



sin f |1)« 



(27) 



This yields a state which is also a ground state of h b z = 
Z® 2 terms at the boundary, but not of the correspond- 
ing X® 2 boundary terms; in other words, charges (Pauli 
Z errors) can condense at the boundaries of the cylin- 
der [27]. The full Hamiltonian — including the h b z terms 
at the boundary — has a two-fold degenerate ground state 
which is topologically protected, and where the logical X 
and Z operators are a loop of Pauli X's around the cylin- 
der and a string of Pauli Z's between its two ends (where 
they condense), respectively. 

To compute pi, we start by considering the PEPS on 
the cylinder without the boundary conditions (27), i.e., 



with open virtual indices at both ends (labelled B and 
B r ). Cutting the cylinder in the middle leaves us with 
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with the joint reduced density operator for the vir- 
tual spaces at the boundary, B (or B'), and the cut, L 
(or R). From (26), one can readily infer that the transfer 
operator for a single tensor is I® 4 + X® 4 , and thus, 



1 Q 



X 



>X* 



(28) 



where the two tensor factors correspond to the B (B 1 ) 
and L (R) boundary, respectively. Imposing the bound- 
ary condition \xe)(xo\i Eq. (27), at B (B'), we find that 
(up to normalization) 

Pi oc (1 + sin 2 0) t® N " + (2 sin 0) X^ , 

which is the thermal state p£ oc exp[— /3 e ffiJg] of Hi = 
—sign (sin 0) X® Nv at an effective inverse temperature 



/?eff 





2sin<9 




tanh 1 








.1 + sin 2 (9_ 





The fact that ijg acts globally is a signature of the 
topological order, and comes from the fact that mea- 
suring an X loop operator gives a non-trivial outcome 
(namely sin#). Note that the entropy S(pe) increases 
by one as l//3 e ff g° es from zero to infinity. This can 
be understood as creating an entangled pair of charges 
| vac) + /(/3 e ff) l c 7 c *) across the cut, thereby additionally 
entangling the two sides by at most an ebit, and subse- 
quently condensing the charges at the boundaries. 

Instead of considering o^, one can also see the topolog- 
ical order by looking at <jbl'- It is the zero-temperature 
state of a completely non-local Hamiltonian X® Nv 
X® Nv which acts simultaneously on both boundaries in a 
maximally non-local way; this relates to the fact that the 
expectation values of any two X loop operators around 
the cylinder are correlated. 

Let us point out that systems with conventional long- 
range order behave quite differently, even though they 
also exhibit correlations between distant boundaries. 
Consider the spin-| Ising model without field, which has 
a PEPS tensor 

^ai,a 2 ,a3,a4 — ^ai,a 2 ^0:2,0:3 ^0:3, 0:4 " 

The resulting local transfer operator is |O)(O| 04 + |l)(l| 04 , 
and thus, 

a BL = 10X01®"' + |1)(1|^ . 

By imposing boundary conditions at £?, one arrives at 

^ = sin#|O)(O| ^+cos#|l)(l|^ , 

which is the thermal state of of the classical Ising Hamil- 
tonian 



log tan 9 

2pN v 



for j3 — >• oo. Thus, for the Ising model, p£ is described 
by a local Ising Hamiltonian, rather than a completely 



non-local interaction as for Kitaev's code state. The 
same holds true for gbl, which is the ground state of 
a classical Ising model without field: while it has correla- 
tions between the two boundaries, they arise from a local 
(i.e., few-body) interaction coupling the two boundaries, 
rather than from terms acting on all sites on both bound- 
aries together. Correspondingly, the long-range correla- 
tions in the Ising model can be already detected by mea- 
suring local observables, instead of topologically nontriv- 
ial loop operators as for Kitaev's code state. 
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FIG. 14: (Color on line) (a) Checkerboard decomposition in 
the Kitaev code. Spin-| are represented by (red) dots at the 
vertices of the square lattice. X and Z operators act on the 
4 spins of each type of (shaded and non-shaded) plaquettes. 
(b) PEPS representation of the Kitaev code (see text). 



VII. CONCLUSIONS AND OUTLOOK 

In this paper, we have introduced a framework which 
allows to associate the bulk of a system with its bound- 
ary in the spirit of the holographic principle. To this 
end, we have employed the framework of Projected En- 
tangled Pair States (PEPS) which provide a natural map- 
ping between the bulk and the boundary, where the latter 
is given by the virtual degrees of freedom of the PEPS. 
This framework allows to map the state of any region to 
a Hamiltonian on its boundary, in such a way that the 
properties of the bulk system, such as entanglement spec- 
trum or correlation length, are reflected in the properties 
of the Hamiltonian. Since our framework also identifies 
observables in the bulk with observables on the boundary, 
it establishes a general holographic principle for quantum 
lattice systems based on PEPS. 

In order to elucidate the connection between the bulk 
system and the boundary Hamiltonian, we have numer- 
ically studied the AKLT model and the Ising PEPS. 
We found that the Hamiltonian is local for systems in 
a gapped phase with local order, whereas a diverging in- 
teraction length of the Hamiltonian is observed when the 
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system approaches a phase transition, and topological 
order is reflected in a Hamiltonian with fully non-local 
interactions; thus, the quantum phase of the bulk can be 
read off the properties of the boundary model. 

Our holographic mapping between the bulk and the 
boundary in the PEPS formalism has further implica- 
tions. In particular, the contraction of PEPS in numeri- 
cal simulations requires to approximate the boundary op- 
erator by one with a smaller bond dimension, which can 
be done efficiently if the boundary describes the ther- 
mal state of a local Hamiltonian, i.e., for non-critical 
systems. Also, since renormalization in the PEPS for- 
malism requires to discard the degrees of freedom in the 
bond space with the least weight [25] , the duality allows 
to understand real space renormalization in the bulk as 
Hamiltonian renormalization on the boundary. 

Our techniques can also be applied to systems in higher 
dimensions, and in fact to arbitrary graphs, to relate the 
boundary of a system with its bulk properties. The map- 
ping applies to arbitrary regions in the lattice, such as 
simply connected (e.g., square) regions used for instance 
for the computation of topological entropies. Also, relat- 
ing the bulk to the boundary using the PEPS descrip- 
tion can be generalized beyond spin systems by consid- 
ering fermionic or anyonic PEPS [28 , as well as conti- 
nous PEPS in the case of field theories [29 , 30 . Finally, 
when studying edge modes, the one-dimensional system 



which describes the physical boundary is given by a Ma- 
trix Product Operator acting on the virtual boundary 
state, and thus, the relation between bulk properties and 
the virtual boundary implies a relation between the prop- 
erties of the bulk and its edge modes physics. 
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